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ABSTRACT 


Understanding  turbulence  degradation  of  electromagnetic 
wave  propagation  is  essential  for  efficient  operation  of  laser 
weapons,  target  designators,  and  imaging  systems.  Random 
atmospheric  refractive  index  inhomogeneities  alter  the  phase 
and  amplitude  of  electromagnetic  waves. 

This  thesis  attempts  to  model  atmospheric  turbulence 
effects  by  using  filtered  Gaussian  phase  screens  to  represent 
the  random  nature  of  refractive  index  changes.  The  simulation 
uses  two-dimensional  512  x  512  fast  Fourier  transform  (FFT) 
techniques  with  extended  Huygens-Fresnel  principles  performed 
on  a  desk  top  computer. 

Simulation  verification  was  accomplished  by  comparing 
calculated  and  theoretical  spatial  coherence  lengths,  po. 
Phase  only  screens  produced  coherence  lengths  that  were  30% 
larger  than  theoretical  values.  By  using  random  phase  and 
amplitude  screens,  the  calculated  coherence  lengths  agreed  to 
within  3%  of  theoretical  values.  Saturation  of  the  normalized 
intensity  variance,  o2/I2,  occurred  for  increasing  turbulence 
using  a  single  phase-amplitude  screen. 
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I.  INTRODUCTION 


Turbulence  in  the  atmosphere  has  a  detrimental  effect  on 
r,  the  propagation  of  electromagnetic  waves  at  optical 

wavelengths.  Minute  changes  in  the  refractive  index  of  the 
atmosphere  along  the  path  of  the  wave  cause  perturbations  in 
the  wavefront  in  both  intensity  and  phase.  Consequently, 
random  atmospheric  density  fluctuations  degrade  the 
performance  of  optical  imaging  systems.  Understanding  this 
phenomenon  is  essential  for  future  endeavors  in  laser  weapons 
and  designators,  as  well  as  imaging  systems. 

This  thesis  attempts  to  model  the  effects  of  atmospheric 
turbulence  using  the  extended  Huygens-Presnel  principle.  This 
theory  considered  both  Fraunhofer  and  Fresnel  propagation 
although  the  simulation  focused  on  a  single  phase  screen  using 
the  Fraunhofer  approximation. 

A  Rytov  approximation  filtered,  Gaussian  distributed, 
random  phase  and  amplitude  screen  simulated  the  effects  of 
turbulence  on  an  optical  beam.  The  magnitude  of  turbulence 
of  this  screen  was  proportional  to  the  index  of  refraction 
structure  parameter,  CBa. 

The  simulation  code  used  two-dimensional  fast  Fourier 
Transform  (FFT)  algorithms  extensively  to  model  optical  wave 

!  * 

propagation  in  the  Fraunhofer  regime. 
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Accuracy  of  the  simulation  was  checked  by  comparison  of 
calculated  and  theoretical  coherence  lengths,  p<> ,  which  were 
obtained  from  the  e-1  point  of  the  atmospheric  mutual 

* 

coherence  function  (MCF) .  Another  check  verified  that 
intensity  variance  saturation  occured  for  an  increase  in  /» 

turbulence,  as  measured  experimentally  and  predicted  by 
theory. 
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II .  BACKGROUND 


A.  STATISTICAL  DESCRIPTION  OF  ATMOSPHERIC  TURBULENCE 
1.  Stationarity ,  Homogeneity,  and  Isotropy 

Most  theorems  involving  random  processes  require  that 
the  functions  satisfy  certain  restrictions.  "Stationarity" 
implies  that  the  mean  value  of  the  random  function  does  not 
change  with  time.  "Homogeneity"  implies  that  translations  in 
x,  y,  and  z  coordinates  (a  Galilean  coordinate  transformation) 
do  not  affect  variables  of  the  random  function.  "Isotropy" 
means  that  the  random  function  is  independent  of  any 
rotational  coordinate  transformation.  [Ref.  1] 

The  last  two  conditions  imply  that  the  statistical 
representation  of  the  random  function  for  two  points  f(ri)  and 
f(rz)  depend  only  on  the  magnitude  of  their  separation 
|  ri  -  rz  j .  [Ref .  1] 

If  these  three  restrictions  only  hold  for  a  local 
volume  of  space,  L3 ,  the  random  field  is  said  to  be  locally 
stationary,  locally  homogeneous,  and  locally  isotropic  for 
|  ri  -  rs  j  <  L. 

Turbulence  of  the  atmosphere  is  a  process  that  can  be 
best  described  by  a  random  function.  Minimal  violations  of 
stationarity,  homogeneity,  and  isotropy  occur  if  the 
neighborhood  is  sufficiently  small,  less  than  some  outer  scale 
length  L.  [Ref.  1] 
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2.  The  Structure  Function 


Because  of  their  nature,  random  functions,  f ,  can  best 
be  described  stochastically.  The  simplest  moment  of  f  is  the 
mean  value,  <f>. 

Another  common  parameter  used  to  represent  a  random 
function  is  the  variance,  o* ,  defined  as 

o*  —  <  {  f (r)  -  <f(r)>  )*  >.  (2.1) 

In  a  physical  random  field,  such  as  atmospheric 
turbulence,  the  mean  value,  <f(r)>,  is  difficult  to  ascertain 
because  of  trends  in  the  mean  as  well  as  variations  in  spatial 
position,  vertically  and  horizontally.  Therefore,  the  mean 
and  variance  are  awkward  parameters  for  describing  a  real , 
random,  atmospheric  field.  To  resolve  this  difficulty, 
Kolmogorov  introduced  the  structure  function 

Df  (r)  =  Df(r*-rt)  =  <[f(rt)  -  f(r2)]2>,  (2.2) 

which  resembles  the  variance  at  first  glance.  [Ref.  2]  The 
structure  function  is  an  ensemble  average  taken  over  all 
possible  point  pairs  ri  and  rj 2 .  [Ref.  2] 

Through  dimensional  analysis,  and  assuming  isotropy, 
Kolmogorov  deduced  a  relation  between  the  structure  function 
and  the  distance  between  the  sampling  points, 

Dt  (r)  «  Cf*  r8/3  (2.3) 
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which  holds  over  a  United  volume,  called  the  inertral 
subrange.  [Re£.  2]  Although  developed  for  the  velocity  field 
v,  Equation  2.3  *lso  holds  for  certain  atmospheric  functions 
called  conservative  passive  additives. 

Conservative  passive  additives  are  quantities  that 
play  no  part  in  the  turbulence  of  the  medium.  Temperature  is 
not  such  a  quantity,  but  potential  temperature,  which  corrects 
for  the  adiabatic  change  in  temperature  with  altitude,  is. 
Therefore,  the  relation 

Dr (r)  =  Ct*  r2/a  (2.4) 

holds  within  the  inertial  subrange  where  T  is  the  potential 
temperature . 

In  examining  problems  of  optical  propagation,  the 
index  of  refraction  is  the  essential  parameter.  To  express 
the  index  of  refraction  as  a  conservative  passive  additive, 
it  is  useful  to  manipulate  the  following  relation 

n  =  1  +  77.6(1  +  7.52  x  10-®/X2)  (P/T)  x  10-«  (2.5) 

where  X  is  the  wavelength  of  the  radiation  in  urn,  P  is  the 
atmospheric  pressure  in  mb,  and  T  is  the  atmospheric 
temperature  in  K.  [Ref.  3]  The  differential  of  Equation  2.5 
is 

dn  *  -  79  x  10-*  P  dT/Ta  (2.6) 
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for  red  light  (X  ■  0.6  pm)  after  ignoring  the  dP  tern 
(isobaric  turbulence) .  Since  potential  temperature  is  a 
passive  additive,  the  index  of  refraction  is  also.  Using 
Equation  2.6  to  express  the  index  of  refraction  structure 
parameter  in  terms  of  the  temperature  structure  parameter 
gives 

Da  (r )  -  Ca*  r*/a  =  (79  P/T2  x  10-*)2  CT2  r2'8  (2.7) 

[Ref.  3]. 

3.  The  Correlation  Punction  and  Power  Spectral  Density 
Expanding  Equation  2.2  gives 

Dr (r)  =  <f2 (ri ) >  -  2  <  f { r i ) f ( r» ) >  +  <f2(r2)>.  (2.8) 

But  if  we  assume  an  homogeneous  medium,  then  <  f2(r)  >  is  the 
same  for  any  point  r,  assuming  <f(r)>  *  0.  The  structure 
function  becomes 

Di  (r)  =  2  Bt  (0)  -  2  Bt  (r),  (2.9) 

where  Tatarski  defined  the  correlation  function,  Bt ,  as 

Bf (r)  -  <  f (ri )  f* (r2 )  >  (2.10) 

in  a  stationary,  homogeneous,  and  isotropic  medium.  [Ref.  2] 
Using  stochastic  Pourier-Stielt jes  integrals,  Tatarski 
developed  a  relation  for  the  three-dimensional  power  spectral 
density  function,  1>(f),  as  the  Fourier  transform  of  the 
correlation  function.  Writing  this  in  freqency  if)  form  as 
opposed  to  radian  freqency  (K)  form  preferred  by  Tatarski, 


6 


exp  (-2nif*r)  B(r)  d8r, 


(2.11) 


and  conversely. 


exp  (-2nif*r)  $(f)  d3f. 


(2.12) 


[Ref.  2] 

Using  the  fact  that  the  correlation  function  is  even,  and 
using  relations  (2.7)  and  (2.9),  the  structure  function  for 
the  Kolmogorov  power  spectral  density  is 


*«  (f)  =  1.303  C»*  f  -*1/3  , 


(2.13) 


for  the  Kolmogorov  power  spectral  density.  This  equation  is 
analogous  to 


(JO  =  0.033  Cn*  JT"* 1  /a  , 


(2.14) 


an  expression  developed  by  Tatarski  in  K  space  that  is  seen 
in  many  publications.  [Ref.  2] 


B.  ELECTROMAGNETIC  PROPAGATION  THROUGH  TURBULENCE 
1.  The  Wave  Equation 

Central  to  any  study  of  electromagnetic  propagation 
are  the  four  equations  of  Maxwell,  presented  here  in  Gaussian 
units 


V  •  H  =  0, 

V  x  H  =  -  ikn2  E , 


(2.15) 

(2.16) 


V 


(n*E) 


0, 


(2.17) 


V  x  E  *  ikH,  (2.18) 

for  the  assumptions  that  the  medium  has  zero  conductivity, 
unit  magnetic  permeability,  and  sinusoidal  dependence. 

Taking  the  curl  of  Equation  2.18  gives 

V  x  (V  x  E)  *  V  x  (ikH).  (2.19) 

Using  a  vector  identity  and  Equation  2.16,  Equation  2.19 
becomes 

-V2 E  +  V(V  •  E)  =  k*n*E.  (2.20) 

Expanding  Equation  2.17  gives  a  relation  for 

V  •  E,  which  when  inserted  into  Equation  2.20  gives  the  wave 
equation 

V2E  +  k2n2E  +  2V[E  •  V(log  n) ]  =  0.  (2.21) 

If  the  wavelength  of  the  electromagnetic  wave  is  small 
compared  to  the  dimensions  of  the  refractive  inhomogeneities, 
then  the  2V[E  •  V(log  n) ]  term  is  negligible.  In  this  case, 
the  wave  equation  reduces  to  the  simple  form 

V2  E  +  k*n2E  =  0.  (2.22) 

[Ref.  3] 

2.  The  Born  Approximation 

One  method  of  solving  the  wave  equation,  employed  by 
both  Tatarski  [Ref.  2]  and  Clifford  [Ref.  3),  is  the  method 
of  small  perturbations  or  the  Born  approximation.  This  method 
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expands  the  electric  field  and  index  of  refraction  as  a  power 
series 

E  =  Eo  +  Ei  +  •••,  (2.23) 

n  *  1  +  m  +  (2.24) 

where  Eo  is  a  constant  in  time,  and  Ei  and  ni  are  time 
varying. 

Inserting  these  two  relations  into  Equation  2.22  and 
equating  same  order  terms  gives 

V2  Eo  +  k2Eo  =  0,  (2.25) 

V2  Ei  +  k2Ei  +  2mk2Eo  =  0,  (2.26) 

ignoring  all  second  and  higher  order  terms. 

Assuming  that  the  electromagnetic  wave  propagates  in 
the  z-direction  and  E  o  *  exp(ikz)  ,  then  Equation  2.26  becomes 

V2Ei  +  k2Ei  +  2ni  k2exp  (  ikz )  =0.  (2.27) 

Clifford  solved  this  wave  equation  for  Ei  by  utilizing 
a  Green’s  function.  The  solution  is  the  convolution  of  a 
plane  wave  Green's  function  with  the  source  term,  shown  here 
in  scalar  form 


Ei  (r)  =  —  / 

4n  J  v 


exp[ik |r-r ' j ] 


[2k2  ni(r’)  exp(ikz')]  d3r, 


(2.28) 


where  r  is  the  position  of  a  specific  point  in  the  image  plane 
and  r’  is  the  postion  of  specific  source  point.  [Ref.  3] 
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Figure  (2.1)  shows  the  coordinate  system  used  to  express  the 
solution  to  the  wave  equation  for  propagation  between  planes 
z'  and  p'  to  z  and  p. 

For  normal  propagation  situations,  |z-z' |  =  R  >>  J />  ~ 
5  'I  so  that  Equation  2.28  expands  into  the  form 

1  r  exp{  ikR  [1  +  \  p  -  p  •  \*/2R‘  +  •••  J| 

Ei  (?)  =  —  /  - 

4n  J  v  R  [1  +  \  p  ~  p  '  j*  /2R*  +  •••] 

X  [2k2  ni(r')  exp(ikz’)]  d2r,  (2.29) 


using  the  binomial  expansion.  Dropping  those  terms  consistent 
with  the  Fresnel  approximation  in  Huygens-Fresnel  optics 
theory,  this  equation  reduces  to 


Ei (r)  « 


k2  exp(ikz) 


2nR 


I 


exp  {  ik  |  p  ~  p  '  | 2  /2R I  m(r')  d^r, 


(2.30) 


which  is  known  as  the  Fresnel  diffraction  formula.  [Ref.  3) 
3.  The  Rytov  Approximation 

Another  method  of  solving  the  wave  equation,  also  used 
by  both  Tatarski  and  Clifford,  is  the  method  of  smooth 
perturbations  or  the  Rytov  approximation.  This  technique 
assumes  that  the  electric  field  is  of  the  form 

E  =  exp(^)  =  exp(X+iS)  *  A  exp(iS).  (2.31) 
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Using  power  series  expansions  for  the  electric  field, 
Equations  2.25  and  2.26  still  hold.  Placing  Equation  2.31 
into  these  two  gives 


(V*E)/E  ♦  k*n*(r)  -  0.  (2.32) 

Tatarski  shows  that  this  method  gives  the  sane  results  as  the 
aethod  of  snail  perturbations  but  the  range  of  validity  is 
larger  than  for  the  Born  approximation.  Equation  2.26  holds 
as  long  as  V*B  is  snail  compared  to  A.  [Ref.  2] 


C.  HUYGENS-FRESNEL  THEORY 

Another  way  of  obtaining  a  relation  for  the  electric  field 
amplitude  at  the  observation  plane,  is  by  extending  the  ideas 
of  the  Huygens-Fresnel  theory  to  include  a  random  phase  term. 
The  Huygens-Fresnel  theory  offers  this  solution  for  the 
electric  field  seen  in  the  observation  plane 


E(r)  *  -  ik  /  E(r')  exp ( ik |  r-r ' | )  Ap  (2.33) 
2n  J s  | r-r '  | 

[Ref.  5] 

Lutomirski  and  Yura  develop  an  extension  of  the  Huygens- 
Fresnel  idea  which  incorporates  a  random  phase  term  which  is 
equivalent  to  the  form  used  in  the  Rytov  approximation, 
exp(^),  where  f  is  a  complex  phase.  [Ref.  4]  This  extended 
integral  is 


E(r ) 


E  (r  ’ )  exp(^(r))  exp  ( ik  (r-r  •  | ) 

jr-r ' | 


12 


(2.34) 


In  the  geometrical  optics  limit,  f  (r)  becomes  the  optical 
phase 


(r)  «  ik 


dz. 


(2.35) 


Expanding  into  a  power  series.  Equation  2.35  reduces  to 
Equation  2.30  as  developed  by  Tatarski  and  Clifford. 


D.  MUTUAL  COHERENCE  FUNCTION 

In  trying  to  get  a  quantitative  measure  for  turbulence, 
Lutomirski  and  Yura  [Ref.  4]  examined  the  average  intensity 
of  the  wave  at  the  observation  point.  The  average  intensity 
at  r  is  defined  as 

<1  (r )  >  *  <  E(ri)  EMra)  >.  (2.36) 


From  Equation  2.34: 
<I(r)>  -  <  (k*/4n* ) 


E(rt')BMrt')  exp[^(?i)  +  ^lr«)] 


X  exp[ik(  jri-ri '  |  -  {ra-ra'|)]  d^Oi  d^a  > 
|ri-?»’J  jrt-rr  ’  | 


(2.37) 


Picking  out  only  the  time-dependent  portion  of  this  result, 
the  entire  integral  reduces  to  evaluating 


<  exp[  (ri)  +  W *  ( r  a )  J  >. 


(2.38) 


This  term  is  the  atmospheric  mutual  coherence  function  (MCF) 
which  contains  the  elements  of  the  propagation  through 
turbulence . 

Lutomirski  and  Yura  [Ref.  4]  relate  the  atmospheric  MCF 
to  the  wave  structure  function,  D,  by  the  relation 
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(2.39) 


MCF(p  )  =  <exp[ (ri )  +  ^  Mr2 )  3  >  *  exp[-D {p  ) / 2 ] 

where  D (p)  *  D x(p)  +  Ds  ( p  )  from  the  Rytov  relation  ^  =  X+iS 
or 

MCF ( p  )  =  exp[-(  p/po)*'3)  ,  (2.40) 


where 


[  1.46  k2 


Cn2  dz 


(2.41) 


and  k  is  the  wavenumber.  Cn2  is  the  refractive  index 
structure  parameter  defined  by  Equation  2.7,  and  R  is  the 
distance  of  propagation  through  the  atmosphere.  The  coherence 
length,  p  o ,  is  defined  as  the  distance  where  transverse 
spatial  coherence  of  the  wave  drops  to  the  e-1  .  [Ref.  1] 

E.  THE  DIFFRACTION  INTEGRAL 

Like  Lutomirski  and  Yura,  Roberts  [Ref.  6]  begins  with  the 
Huygens-Fresnel  integral  for  the  propagation  of  light  waves 
after  making  the  Fresnel  approximation  (|z-z'|  =  R  >>  \p~p' | ) 


E  ( r ) 


-ik 

2nR 


exp  {  ik  [R2  +  |  p  p  |2)‘/2 . 


(2.42) 


In  this  expression,  k  is  the  wavenumber  and  the  notation  is 
the  same  as  shown  in  Figure  (2.1). 

Employing  the  standard  practice  of  expanding  by  the 
binomial  expansion  and  neglecting  smaller  terms  gives 
E(r)  =  -ik  exp[-ikRj  f  dp  ’  E(r’)  exp{ik|  p'~p  |2  I  • 


2nR 


2R 
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(2.43) 


Since  the  exponential  term  outside  the  integral  doesn't  affect 
intensity  measurements,  it  can  be  neglected.  Expanding  the 
quadratic  term  gives  the  relation 


E (r )  =  -ik  /  dp  '  E(r')  exp{ik  [  p  *2-2 p  ' • p +  p 2 ] I 
2nR  /  s  2R 

=  -ik  exp[ik  p'2 1  I  dp'  E(r')  exp [ ik  t]  expf-ikiQ  »  p  '  1  . 
2nR  2Rr  J  s  2R  'r 

(2.44) 

This  is  called  the  convolution  form  of  the  Fresnel  diffraction 
integral  which  differs  from  the  Fraunhofer  approximation  only 
by  the  quadratic  phase  term  exp[ik^j_2].  [Ref.  6] 

If  we  denote  the  Fourier  transform  of  E(r)  as  E(f),  then 


E(  f)  -Id p  exp(-2nif*/3  )  E(r) 


(2.45) 


Inserting  E(r)  as  defined  in  Equation  2.43,  gives 

E(  f)  =  -ik  dp  '  E  (r  '  )  I  dp  exp  ( -2n if-p  )  exp  I  ik  |  p  '  ~P  I 
2nR  J  s  J  s  2R 


(2.46) 


Changing  variables  to  q  =  p  '  -  p  ,  the  electric  field 

spectrum  becomes 


E(f)  =  -ik  f  dp  '  E(r')  exp(-2n if-5 
2nR  Js  r 


) 


/  dq  exp(-2nif*q)  exp(ik  q2 ) 
Is  2R 


(2.47) 


The  q  integral  is  shown  to  be  the  Fourier  transform  of  the 
Gaussian  function 
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2niR  exp  ( ~2n2  iRf*  ) 
k  k 


(2.48) 


[Ref.  6]  Substituting  this  function  into  Equation  2.47  gives 

E(r')  exp(-2n if'jo),  (2.49) 


E(f)  *  exp  (-2n2  iRf* ) 
k 


which  is  the  inverse  transform  of  what  we  sought. 


E  (r ) 


■/ 


df  exp(2ni/*r)  exp (~2n8-iRfa  ) 

k 


f«P  ’ 


E (r ' )  exp (-2nif • p  ) 


(2.50) 


This  is  the  transfer  form  of  the  Fresnel  diffraction  integral. 
[Ref.  6] 

Equations  2.44  and  2.50  are  two  different  forms  of  the 
diffraction  integral  that  are  useful  for  two  types  of  wave 
propagation.  The  convolution  form,  Equation  2.44,  is  best 
suited  for  long  distance  propagation  (Fraunhofer  regime)  since 
R  is  in  the  denominator  of  the  exponential  term.  The  transfer 
form.  Equation  2.50,  is  best  suited  for  short  distance 
propagation  since  R  is  present  in  the  numerator  of  the 
exponential  term.  [Ref.  6] 


16 


III.  NUMERICAL  SIMULATION 


A.  SIMULATION  HARDWARE  AMD  SOFTWARE  SPECIFICS 

This  numerical  simulation  modeled  the  propagation  of  an 
electromagnetic  wave  through  a  random  turbulent  medium.  The 
procedures  for  this  simulation  required  the  creation  and 
manipulation  of  multiple  two-dimensional,  512  by  512,  single 
precision,  phase  screens,  requiring  a  minimum  memory  of 
several  megabytes . 

A  COMPAQ  Deskpro  80386-20  computer  configured  with  a  16 
megabyte  RAM  and  64  megabyte  hard  drive  met  hardware 
requirements.  The  computer  also  had  EGA-VGA  graphics  and  both 
HP  Laser  Jet  II  and  Panasonic  KX-P1092i  printers.  A  Weitek 
1167  math  coprocessor  increased  execution  speeds  of  the  20  Mhz 
80387  coprocessor  by  a  factor  of  3-4. 

The  software  requirements  for  writing  the  simulation  code 
were  met  with  the  MicroWay  32  bit  NDP  FORTRAN-386  compiler 
and  the  Phar  Lap  operating  system  extender.  It  was  chosen 
over  the  Science  Applications  International  Corporation's  SVS 
Fortran-386  compiler  which  had  10-20%  faster  program  execution 
times  but  numerous  bugs  present  in  the  graphics  routines. 

Unfortunately,  due  to  limitations  of  version  1.4e  of  the 
NDP  Fortran-386  compiler,  the  full  graphic  capabilities  of  VGA 
were  not  supported,  so  EGA  graphics  were  used  for  screen 
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displays.  Hardcopy  graphs  and  plots  were  accomplished  using 
Grapher  and  Surfer  programs  from  Golden  Software. 

B.  REQUIRED  ALGORITHMS 

For  the  successful  operation  of  such  a  model  both  the 
random  number  generator  and  the  two-dimensional  fast  Fourier 
transform  (FFT)  routines  needed  speed  and  accuracy.  These 
will  be  discussed  in  detail  in  the  following  sections. 

1.  Random  number  generator 

The  effects  of  turbulence  on  the  propagation  of  an 
electromagnetic  wave  are  of  a  random  nature.  To  emulate  this 
random  nature  required  a  pseudo-random  number  generator 
algorithm.  The  minimum  criteria  for  a  quality  random  number 
generator  was  that  the  products  be:  (1)  sufficiently  random, 
(2)  reproducible,  and  (3)  rapxdiy  obtained. 

The  NDP  FORTkAN  compiler  did  not  come  equipped  with  a 
random  number  generator,  so  an  algorithm  of  the  linear 
congruent  type  especially  designed  for  32  bit  machines  was 
implemented  as  a  subroutine  attached  to  the  main  program. 
This  had  the  benefit  of  not  binding  this  portion  of  the  code 
to  a  specific  machine  setup.  This  random  number  generator 
algorithm  came  from  Dr.  Milne,  who  obtained  it  from  class 
notes  prepared  by  Dr.  Harrison  [Ref.  7],  This  algorithm 
appears  as  subroutine  RAN  (Appendix  A)  in  the  simulation  code. 

The  first  criteria,  that  the  generator  produce  numbers 
sufficiently  random,  was  hard  to  quantify.  A  series  of  tests 
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exercised  the  statistical  independence  of  succeeding  numbers. 
There  are  many  such  tests  that  exist,  the  most  common  of  which 
is  the  "Chi-Square"  test.  The  following  is  a  list  of  all 
tests  used  to  check  the  quality  of  the  random  number  generator 
algorithm: 

1.  Frequency  Test.  This  tests  the  uniformity  of  the 
sequence  of  numbers . 

2.  Serial  Test.  This  tests  the  two-dimensional  uniformity 
of  the  sequence  of  numbers . 

3.  Runs  Up  and  Down  Test.  This  tests  to  see  how  long 
sequences  of  random  numbers  are  that  either  go 
successively  up  or  down. 

4.  Laqqed-Product  Test.  This  tests  for  correlations 
between  random  numbers  Ri  and  Ri*j  where  j  is  the  given 
lag  value.  A  lag  of  3  is  especially  difficult  to  pass. 

5.  Repeat  Test.  This  simple  qualitative  check  is  to 
determine  how  long  the  sequence  of  random  numbers 
becomes  before  the  first  number  is  repeated. 

6.  Scatter  Plot  Test.  This  is  another  qualitative  test 
where  random  numbers  plotted  on  a  x-y  plane  are  visually 
checked  for  correlation. 


These  six  tests  were  performed  for  a  variety  of  seed 
values  to  determine  which  gave  the  best  performance.  Results 
of  the  first  five  tests  for  two  good  seeds  can  be  found  in 
Table  3.1. 

Choosing  the  same  initial  seed  value  gave 
reproducibility  of  the  random  numbers.  The  algorithm's  simple 
two  lines  of  code  assured  rapid  execution. 
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TABLE  3.1  RANDOM  NUMBER  GENERATOR  TEST  DATA 


This  table  provides  the  results  of  statistical  tests  of 
the  random  number  generator  algorithm.  The  Xs  values  came 
from  the  Chi-Square  table  in  Bevington  [Ref.  11] 

The  results  of  the  Lagged  Product  test  correspond  to  the 
theoretical  mean,  ut  ,  the  calculated  mean,  \i,  the  theoretical 
standard  deviation,  ax,  and  the  calculated  standard  deviation. 


Test 

Seed 

Degrees  of 
Freedom 

X* 

Probability 

45739 

9 

■DEI 

1 

91 

% 

Frequency 

MM 

■ 

94377 

9 

■ 

1 

82 

% 

45739 

20 

16.6 

70  %  , 

Serial 

94377 

20 

10.2 

96 

% 

Runs 

45739 

6 

■WEI 

■ 

2 

% 

Up  and 

Mi 

Down 

94377 

5 

■ 

80 

% 

Test 

Seed 

P* 

P 

or 

o 

45’739 

0.250 

0.230 

0 

.088 

0.092 

Lagged 

Product 

94377 

0.250 

0.265 

0 

.088 

0.086 

2.  Fast  Fourier  Transform 

Since  the  most  often  used  algorithm  in  the  simulation 
was  the  Fast  Fourier  Transform  (FFT) ,  efficiency  was  critical. 
Two  different  routines  were  compared. 

The  first  FFT  routine  considered  was  a  package  offered 
by  NDP  of  one-  and  two-dimensional  machine  language  FFT 
routines  designed  for  easy  linking  by  the  NDP  FORTRAN 
compiler.  This  package  had  routines  available  for  real  and 
complex  arrays. 
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The  second  FFT  routine  considered  was  a  simple  real 
array  subroutine  coded  by  Dr.  Walters  from  a  BASIC 


demonstration  package  provided  by  Infotek. 

Each  FFT  routine  was  compared  for  numerical  output  and 
speed  of  operation  for  the  one-dimensional  case.  As  expected, 
all  routines  were  virtually  identical  in  numerical  output. 
However,  speed  of  operation  differed  widely,  as  seen  in  Table 
3.2. 


TABLE  3.2  FFT  SPEED  OF  OPERATION  COMPARISON 

This  table  contains  results  of  FFT  execution  time 
comparison  between  NDP  machine  language  routines  and  Dr. 
Walters'  subroutine.  Since  NDP  FORTRAN  does  not  have  a  timing 
function  with  sufficient  accuracy,  execution  times  listed  are 
mean  times  found  from  several  runs,  timed  by  stopwatch. 


Array 

NDP  Routines 

Walters 

Routine 

Size 

CFFT 

RFFT 

w/o  & 

w/  Weitek 

1-D  field  size 

time (s) 

time (s ) 

time (s) 

time (s ) 

2»3 

2.60 

1.45 

2.84 

1.20 

2‘« 

5.27 

2.88 

5.86 

2.32 

213 

10.94 

5.84 

12.27 

4.94 

216 

22.78 

12.05 

25.90 

10.26 

217 

47.61 

25.02 

54.53 

21.55 

The  differences  in  run  time  for  the  two  NDP  machine 
language  routines  were  due  to  the  CFFT  routine  being  more 
cumbersome  in  converting  repeatingly  between  real  and  complex 
arrays.  The  compiled  subroutine  was  significantly  faster  than 
both  CFFT  and  RFFT  NDP  machine  language  routines  unless  the 
Weitek  was  deactivated.  Using  a  subroutine  enclosed  within 


the  main  program  had  the  additional  benefit  of  not  being 
machine  dependent.  Therefore,  it  was  decided  that  Dr. 
Walters'  subroutine  would  be  used  for  all  PFTs  needed  in  the 
simulation  code.  Conversion  of  this  subroutine  for  two- 
dimensional  use  was  relatively  easy.  The  NDP  two-dimensional 
PFT  routines  never  worked  correctly. 

A  PFT  "breaks  down"  a  function  into  its  cosine  and 
sine  constituents  (real  and  imaginary  terms,  respectively). 
The  range  of  spatial  frequencies  represented  in  a  FFT  depend 
on  the  array  size  chosen.  The  minimum  frequency  is  fain  =  1/W 
and  the  maximum  frequency  is  fBax  =  n/W,  where  W  is  the 
physical  width  of  the  array  and  n  is  the  number  of  sampling 
points  across  the  array.  When  energy  at  frequencies  near 
these  cutoff  values  appeared,  problems  developed  due  to  the 
discrete  nature  of  the  FFT. 

One  problem  exists  in  representing  spatial  frequencies 
smaller  than  fain  due  to  the  high  amplitude,  low  frequency 
"tilt"  of  the  f  filtered  wavefront.  Since  a  tilted 
straight  line  approximates  a  cosine  or  sine  function  for  only 
a  small  region,  it  helps  to  make  the  wavefront  smaller  than 
the  FFT  array  size.  Hence,  to  help  represent  lower 
frequencies  in  the  simulation  coding  accurately,  an  aperture 
array  was  a  user  chosen  variable  smaller  than  the  FFT  array 
size . 

Another  problem  that  manifests  itself  is  "aliasing", 
which  occurs  when  high  frequency  features  of  a  function  are 
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missed  by  having  too  few  sample  points.  This  means  that 
important  information  for  the  function  present  in  the 
frequencies  greater  than  fMax  are  folded  back  and  appear  at 
lower  frequencies.  Choosing  the  distance  between  sample 
points  smaller  than  the  coherence  length  of  the  propagating 
electromagnetic  wave,  found  theoretically  by  Equation  2.39, 
alleviated  this  problem. 


C.  SIMULATION  DESCRIPTION 

1.  User  defined  quantities 

For  operation  of  the  simulation  program,  the  user  had 
the  option  of  varying  several  quantities  directly  from  the 
keyboard: 

1.  Array  size.  This  can  be  varied  up  to  a  maximum  of  512 
by  512  points. 

2.  Aperture  sampling  points.  This  can  be  varied  up  to  a 
maximum  of  512  by  512  points.  Ideally  it  should  be 
smaller  than  the  array  size  to  avoid  "tilt"  effects,  but 
large  enough  to  avoid  "aliasing"  effects.  A  size  one- 
fourth  to  one-half  the  array  size  works  best. 

3.  Aperture  shape .  A  choice  of  square  or  circular  aperture 
was  provided. 

4.  Seed .  The  pseudo-random  number  generator  required  a 
five  digit  input  seed.  Table  3.1  lists  two  quality 
seeds.  All  data  was  obtained  using  the  seed  94377. 

5.  Turbulence  parameter.  Any  value  of  Cnz  above  zero  can 
be  chosen.  Typical  values  of  C»2  were  10"17  through 

10-13  m-2/3. 

6.  Wavelength .  Any  value  of  wavelength  for  the  electro¬ 
magnetic  wave  can  be  chosen.  A  value  of  0.6  pm  was  used 
throughout . 

7.  Distance .  This  was  the  distance  from  the  aperture  to 
the  observer.  This  simulation  used  1  km. 


8.  Aperture  width.  Adjustable  to  ensure  that  "aliasing" 
effects  do  not  occur  for  a  specific  choice  of  C»*. 
Typically,  a  100  nun  aperture  was  chosen  with  a  path 
length  of  1  km. 

For  the  accumulation  of  multiple  data  points  for 
statistical  purposes,  a  slightly  different  version  of  the 
program  read  user  choices  from  a  disk  file. 

2.  Aperture  Mutual  Coherence  Function 

Computing  far-field  diffraction  patterns  of  intensity 
for  an  electromagnetic  wave  through  the  aperture  with  no 
turbulence  gave  the  Mutual  Coherence  Function  (MCF)  of  the 
aperture.  Since  we  assumed  a  point  source  at  infinity,  the 
electric  field  wave  front  was  planar  at  the  aperture.  The 
electric  field  that  passes  through  had  a  value  of  one  within 
the  boundaries  of  the  aperture,  and  a  value  of  zero  outside 
the  boundaries.  The  FFT  subroutine  computed  the  diffraction 
pattern,  the  electric  field  present  at  the  image  plane.  The 
conjugate  square  of  this  electric  field  was  the  intensity 
field  that  would  be  seen  by  an  observer.  Taking  the  inverse 
FFT  of  the  intensity  field  gave  the  aperture  MCF,  which  was 
identical  to  the  autocorrelation  of  the  aperture  function. 

3.  Phase  screen  generation 

Simulation  of  effects  of  turbulence  requires  that  the 
aperture  planar  electric  field  be  multiplied  by  the  phase 
screen.  Two  methods  exist  to  create  the  phase  screens. 
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a.  Classical  method 

In  this  method  n2  Gaussian  distributed  random 
numbers  filled  a  n  by  n  real  array.  Since  the  pseudo-random 
number  generator  RAN  produced  uniformly  distributed  numbers, 
a  subroutine  GAUSS  transformed  them  to  a  Gaussian  distribution 
with  unit  variance,  Knuth  [Ref.  8].  This  subroutine  appears 
in  Appendix  B.  Taking  a  direct  FFT  creates  two  Gaussian 
distributed  arrays,  phaser  and  phasei,  which  are  the  real  and 
imaginary  spectral  amplitudes  of  the  original  random  array. 
Since  these  two  arrays  are  in  spatial  frequency  space, 
multiplication  by  Equation  2.13  occurs.  An  inverse  FFT 
returns  the  filter  to  real  space.  Because  of  the  symmetry  of 
the  filter,  discussed  later,  only  one  array  in  real  space 
results,  containing  all  n2  pieces  of  information  produced  by 
the  subroutine  RAN. 

b.  Martin  and  Flatt6  method 

In  this  second  method  suggested  by  Martin  and 
Flatt§  [Ref.  9],  two  n  by  n  arrays,  called  phaser  and  phasei, 
are  each  filled  with  n2  Gaussian  distributed  random  numbers. 
These  are  considered  the  real  and  imaginary  components  in 
spatial  frequency  space.  Filtering  occurs  as  5*,,/3  then  an 
inverse  FFT  returns  the  array  to  real  space.  Because  no 
information  is  lost  in  filtering,  two  arrays  result,  each  with 
n2  pieces  of  information.  Only  one  of  these  arrays  need  be 
used  for  the  following  simulation  code.  This  method  requires 
twice  the  random  numbers  as  the  previous  method  but  GAUSS 
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produces  two  random  numbers  and  half  as  many  two-dimensional 
FFT  steps  are  required.  For  use  in  Fresnel  propagation,  where 
multiple  phase  screens  exist,  this  method  is  certainly  the 
most  efficient. 

4 .  Filtering 

Phase  screens  need  to  be  filtered  in  frequency  space 
to  ensure  the  proper  power  law  form.  From  Tatarski  [Ref.  1] 
and  Martin  and  Flatty  [Ref.  9],  the  correct  form  for  the 
filtering  function,  4>f  ,  is 

=  2n  k2  5z  $n,  (3.1) 

where  k  is  the  wavenumber  of  the  electromagnetic  wave,  6z  is 
the  slab  thickness,  and  4>n  is  the  power  spectral  density 
function.  Using  Equation  2.13  for  the  power  spectral  density 
function  and  the  definition  of  wavenumber,  k  =  2n/A,  gives 

4>f  =  1.303  ( 2n) 3  (1/A)2  Cn2  5Z  f  -*»'*  (3.2) 

which  is  the  filtering  function  used  in  the  subroutine  FILTER 
in  the  simulation  code.  Arrays  Phaser  and  phasei  are  both 
multiplied  by  the  square  root  of  Equation  3.2  to  obtain  a 
factor  for  the  real  and  imaginary  amplitudes  so  that  the 
intensity  has  the  proper  f  3  power  law. 

Although  this  filtering  scheme  looks  simple,  a  few 
subtle  points  make  the  process  much  more  complicated  than  is 
immediately  evident.  In  fact,  perfecting  this  filtering 
process  constituted  the  major  effort  of  this  thesis. 
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First  of  all,  most  of  the  relations  developed  in  the 
background  references  were  in  terms  of  the  radian  frequency 
K.  The  FFT  routine  used  was  in  terms  of  the  spatial  frequency 
f  =  K/2n.  Consequently,  all  relevant  equations  had  to  be 
expressed  in  this  latter  form,  requiring  the  tracking  down  of 
many  2n  terms. 

Secondly,  due  to  the  two-dimensional  isotropy  of 
turbulence,  the  correct  filtering  scheme  is  of  a  circular 
nature.  This  means  that  the  frequency  seen  in  Equation  3.2 
is  really  f  =  ( f*2  +  fy2)l/2,  which  is  radial  everywhere  but  at 
the  zero  frequency  (DC  term)  where  it  is  has  a  value  of  zero 
[Ref.  9]. 

For  convenience  of  coding,  the  zero  frequency  should 
be  at  the  center  of  the  arrays  phaser  and  phasei.  But  in  the 
FFT  algorithm,  spatial  frequencies  are  arranged  in  a  manner 
that  is  quite  different.  The  zero  frequency  term  is  at  the 
upper  left  corner  of  the  array.  The  Nyquist  "folding" 
frequency  is  present  as  a  cross  that  divides  the  arrays  into 
four  sections.  So,  prior  to  filtering,  a  shuffling  of  data 
points  is  required  as  discussed  by  Brigham  [Ref.  10)  and  as 
seen  in  Figure  3.1.  After  filtering,  an  inverse  shuffle 
returns  all  frequencies  to  their  previous  locations. 

Another  important  consideration  in  the  filtering 
process  was  that  frequency  units  be  correct  according  to  the 
physical  dimensions  of  the  aperture.  This  required  that  the 
frequency  array  have  a  normalization  factor  of  1/NA  where  N 
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is  the  total  number  of  points  across  the  array,  and  A  is  the 
sampling  interval  (physical  units  between  the  aperture  array 
points) . 

5.  Propagation 

The  Huygens-Fresnel  techniques  formed  the  basis  for 
the  simulation's  treatment  of  optical  propagation.  For 
simplicity  this  simulation  only  considered  Fraunhofer 
propagation.  The  addition  of  subroutines  to  handle  the 
dynamics  of  Fresnel  propagation  should  not  be  too  difficult 
other  than  the  need  for  dynamic  scaling  of  the  numerical  mesh 
to  account  for  diffraction  spreading. 

Martin  and  Flatty  provided  a  procedure  to  simulate  the 
propagation  of  a  wave  in  a  turbulent  medium.  [Ref.  9] 


1.  Add  the  random  turbulence  effects  to  the  electric  field 
distribution  by  meshing  the  aperture  electric  field 
function  with  a  random  phase  term  exp(^): 

E  (r '  )  -  E  (r  '  )  exp  (  V  )  .  (3.3) 

2.  Take  the  Fourier  transform  of  the  result  of  Equation  3.3 
to  obtain 

E(f'  )  =  DFT  (E  (r '  )  I  .  (3.4) 

3.  Add  a  term  to  the  result  of  Equation  3.4  to  obtain 

E{  f)  =  E(f')  exp  (  -2niRf2  i  (3.5) 

k 

which  is  the  exact  same  relationship  as  in  Equation 
2.48  that  was  described  by  Roberts  as  the 
convolution  form  of  the  Fresnel  diffraction  integral. 
[Ref.  6] 

4.  Take  the  inverse  Fourier  transform  of  Equation 
3.5  to  get  an  expression  for  the  electric  field 
distribution  at  the  image  plane: 
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E(r)  =  IFT I  E{  f)  ) 


(3.6) 


If  one  were  taking  a  multiple  screen  approach  to 
solving  for  the  propagation  of  the  wave  as  suggested  by  Martin 
and  Flatt6,  the  entire  process  would  be  reiterated  for  as  many 
segments  of  length  R  as  required.  In  this  numerical 
simulation,  only  one  segment  was  used  for  all  propagations. 
Adding  this  iterative  section  to  the  code  would  be  a 
straightforward  modification  if  dynamic  array  sizes  were 
incorporated.  [Ref.  6] 

For  distances  of  propagation  where  the  phase  front 
cannot  be  assumed  planar,  an  additional  quadratic  term  would 
have  to  be  added  at  step  3  for  the  transfer  form  of  the 
Fresnel  diffraction  integral,  as  seen  in  Equation  2.50. 

6.  Atmospheric  Mutual  Coherence  Function 

We  computed  the  atmospheric  mutual  coherence  function 
after  meshing  the  random  phase  screen  with  the  aperture 
electric  field  function  (after  step  1  of  the  Martin  and  Flatt£ 
method)  .  As  with  the  aperture  MCF,  the  next  step  was 
computing  the  diffraction  pattern  using  the  FFT  routine,  which 
represents  the  complex  electric  field  that  an  observer  would 
see  because  of  the  aperture  and  the  turbulence.  The  effects 
of  turbulence  tends  to  spread  out  the  diffraction  pattern  so 
that  resolution  would  be  lost  for  any  optical  imaging  system. 
Taking  the  conjugate  square  of  this  diffraction  pattern  gives 
the  intensity  field,  which  negates  the  additional  exponential 
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term  present  in  Equation  3.5.  Taking  the  inverse  FFT  of  the 
intensity  field  gave  the  total  mutual  coherence  function. 
Dividing  the  total  MCF  by  the  aperture  MCF  leaves  the 
atmospheric  MCF. 

7 .  Coherence  Length 

As  stated  previously,  the  coherence  length,  p*  ,  was 
the  distance  transverse  to  the  direction  of  propagation  where 
the  spatial  coherence  of  the  <E*E>  wave  dropped  to  e_1  . 

For  a  continuous  plot  of  the  atmospheric  MCF,  the 
e-1  distance  could  be  picked  off  the  curve  easily.  Due  to 
the  discrete  nature  of  the  simulation,  the  MCF  curve  was  not 
continuous,  but  a  series  of  points  equal  to  the  choice  of  the 
pixel  width  of  the  aperture  array.  Therefore  this  simulation 
used  a  linear  interpolation  routine  to  pick  out  values  for 
the  coherence  length  as  seen  in  the  results  of  the  following 
section . 


In  order  to  verify  that  this  numerical  simulation  produced 
correct  results,  we  had  to  find  ways  to  check  the  validity  of 
the  code.  In  the  code,  presented  in  Appendix  C,  there  are 
several  checkpoints  where  the  accuracy  of  the  simulation  could 
be  verified.  The  following  sections  discuss  these  checkpoints 
in  increasing  order  of  sophistication. 

A.  APERTURE  HCF 

An  initial  point  to  check  the  accuracy  of  the  simulation 
code  was  to  look  at  the  mutual  coherence  function  (MCE)  of 
the  aperture.  Since  the  code  offered  two  choices  of  aperture 
shape,  square  and  circular,  both  were  verified. 

Using  the  two-dimensional  FFT  routine,  the  diffraction 
patterns  for  100  mm  wide  square  and  circular  apertures  were 
calculated.  The  square  aperture  diffraction  pattern  in  Figure 
4.1a  displays  the  familiar  [(sin  x)/x]*  form  discussed  by 
Hecht  with  the  correct  maximum  value  and  node  spacing.  [Ref. 
5]  The  circular  aperture  diffraction  pattern  in  Figure  4.1b 
displays  the  Airy  pattern  derived  from  the  Ji  Bessel  function 
solution  carried  out  by  Hecht.  [Ref.  5] 

An  inverse  2-D  FFT  of  each  100  mm  wide  aperture  shape 
produced  the  aperture  MCF  (identical  to  aperture 
autocorrelation) .  The  square  aperture  autocorrelation  in 
Figure  4.2a  is  the  predicted  triangle  function.  The 
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Figure  4.1  Irradiance  Diffraction  Pattern  Resulting  From 
A  (a)  Square  Aperture  And  A  (b)  Circular  Aperture. 
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Figure  4.2  Autocorrelation  (Aperture  Mutual  Coherence 
Function)  Of  A  100  mm  Wide  (a)  Square  Aperture  And  A 
100  mm  Wide  (b)  Circular  Aperture 


Length  (meters) 


circular  aperture  autocorrelation  in  Figure  4.2b  ia  aa 
expected  for  a  "top  hat"  function. 

The  algorithm  gave  the  correct  MCP  results  which  verified 
the  simulation  up  through  the  aperture  MCP  point. 

B.  FILTERING 

The  most  crucial  portion  of  the  entire  simulation  code  was 
the  correct  frequency  filtering  algorithm.  A  succinct  way  of 
verifying  this  came  from  examining  the  result  of  taking  the 
logarithm  of  both  sides  of  Equation  3.3 

log{*f|  «  logll . 303 (2n)a  (l/X)aCna5* .f  -»»/«} 
which  is 

log{*f|  -  log{1.303(2n)3  <l/X)*Cn*5zl  -  (11/3)  log!  f)  ,  (4.1) 
an  equation  for  a  straight  line  with  slope  of  -11/3. 

Figure  4.3  shows  a  plot  of  filtering  function  vs. 
frequency  f  with  a  value  10- 13  for  Cna  and  a  wavelength  of  600 
nm  (red  light)  using  the  Martin  and  Flatt6  method  of  filtering 
discussed  earlier.  A  linear  fit  to  the  data  points  down  to 
the  Nyquist  frequency  had  a  slope  of  -3.88.  This  was  a 
deviation  from  the  predicted  value  of  -11/3  (-3.83)  by  1.3%. 
Therefore,  it  appears  that  the  filtering  is  being  accomplished 
correctly.  The  upturn  of  points  occurs  at  the  Nyquist  folding 
point  that  is  present  in  the  array  scheme  for  the  PFT  routine. 

C .  ATMOSPHERIC  MCF 

An  important  place  to  check  results  of  the  simulation  was 
at  the  calculation  of  the  atmospheric  MCF.  Following  the 


35 


Log  of  Phase 


Figure  4.3  Slice  Through  Phase  Screen  Function 
Filtered  By  The  Martin  And  Flatt6  Method  Showing 
Its  Dependence  On  Spatial  Frequency. 

Slope  Of  -11/3  Comes  Directly  From  The  Kolmogorov 
Power  Law.  The  Upturn  At  Higher  Frequencies  Is  Due 
To  Nyquist  Folding  Present  In  The  Phase  Screen  Array  Code. 
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procedure  for  calculating  the  atmospheric  MCF  outlined  in 
section  III,  first  the  diffraction  pattern  was  calculated. 
Figure  4.4a  shows  the  diffraction  pattern  for  a  100  mm  wide 
square  aperture  that  results  for  red  light  (X  «  0.6  vun)  with 
Cn*  ■  10”1*.  Notice  that  the  diffraction  is  only  slightly 
affected  from  what  is  shown  in  Figure  4.1a.  As  the  turbulence 
structure  parameter  Cn*  increased,  the  diffraction  pattern 
began  to  spread  out  more  until  it  no  longer  resembles  the 
[(sin  x)/x]2  form  it  had  without  turbulence,  as  seen  in 
Figures  4.4b,  4.4c,  and  4.4d  for  values  of  Cn*  *  10-18,  Cn*  * 
10-14,  and  Cn*  *  10-la  respectively.  Figure  4.4c  shows  that 
modest  amounts  of  turbulence  displace  the  image  centroid,  due 
to  low  frequency  tilting  of  the  electric  field.  Since  the 
diffraction  pattern  is  analogous  to  the  intensity  that  would 
be  seen  at  the  image  plane,  this  clearly  shows  the  blurring 
of  an  image  commonly  seen  looking  through  a  turbulent 
atmosphere  as  well  as  the  "dancing"  image  effect  at  lower 
levels  of  turbulence. 

The  atmospheric  MCF  was  calculated  for  different  values 
of  turbulence  by  the  method  described  in  section  III  using  the 
classical  method  of  filtering.  These  are  seen  in  Figures  4.5 
for  Cn*  values  of  0,  10"  10-18,  10-14,  and  10-13.  With  no 
turbulence,  the  atmospheric  MCF  is  a  value  of  1.0  everywhere, 
seen  as  a  straight  line.  With  increasing  turbulence,  the 
atmospheric  MCF  curves  fall  off  rapidly  to  zero. 
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(c) 


(d) 


Figure  4.4  Diffraction  Patterns  For  A  100  mm  Wide  Square 
Aperture  For  Four  Different  Values  Of  Turbulence: 

(a)  C.*«10-‘*  (b)  C.*«10-»®  (c)  C«*=10-‘«  (d)  C.*  =  10-»3. 
Displacement  Of  The  Beam  Centroid  Is  Obvious  In  (c)  Due 
To  Low  Frequency  Tilting  Of  The  Electric  Field. 
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D.  COHERENCE  LENGTH 

The  e_l  point  of  atmospheric  MCF  curves  give  the  coherence 
length,  p o.  In  the  simulation,  we  achieved  this  by  a  method 
of  linear  interpolation  as  discussed  previously.  One  should 
not  put  too  much  faith  in  a  statistical  sample  of  one, 
however.  For  this  reason,  a  modification  of  the  simulation 
was  made  to  run  several  samples  of  each  different  set  of  input 
parameters  so  that  a  better  statistical  representation  of 
results  occurred. 

Theoretical  coherence  lengths  were  calculated  from 
Equation  2.41  [Ref.  1]  for  three  different  values  of  Cn2  using 
R=1000  m  and  A“0.6  pm.  Several  runs  of  data  were  produced  to 
investigate  the  effect  of  different  sizes  of  both  the  electric 
field  array  and  the  aperture  array.  Calculated  results  are 
discussed  and  displayed  in  following  pages. 

Figure  4.6  shows  the  trend  for  coherence  length  values  as 
a  function  of  aperture  array  sample  points  using  the  classical 
method  of  filtering  for  Cn2  =  10-18.  Ten  samples  were  obtained 
for  each  data  point.  As  discussed  in  section  III,  aliasing 
makes  the  values  of  p o  too  high  if  there  are  too  few  sample 
points.  Using  aperture  arrays  no  larger  than  one-half  the 
width  of  the  main  array  was  done  to  avoid  tilt  problems.  The 
curves  asymptotically  approach  coherence  length  values  of  4.18 
mm  for  the  512  x  512  array  and  4.49  mm  for  the  256  x  256 
array.  These  represent  errors  of  39.3%  and  49.7% 
respectively . 


40 


I 


S  15 

a 


be 

a 

<D 

V 
O 

ti 

V 
u 

V 

S3 

o 

u 


10H 


0 


1 — i — i — r 


0 


•••••512  by  512  array 
ooooo  256  by  256  array 
- Theoretical  coherence  length 


t — i — i — i — i — | — i — i — i — i — i — i — i — i — i — | — i — n — i — i — p — i — i — i — | 

100  200  300 

Aperture  Sample  Points  (pixels) 


Figure  4.6  Calculated  Coherence  Length  Values  For 
Cn*  *  10“‘3  Using  Different  Aperture  Array  Sizes 
And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 


41 


Figure  4.7  shows  coherence  length  trends  for  Cn*  « 

10-14  also  using  the  classical  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  15.78 
mm  for  the  512  x  512  array  and  18.43  mm  for  the  256  x  256 
array,  representing  errors  of  32.1%  and  54.2%  respectively. 

Figure  4.8  shows  coherence  length  trends  for  Cn*  = 

10-19  also  using  the  classical  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  61.99 
mm  for  the  512  x  512  array  and  76.27  mm  for  the  256  x  256 
array,  representing  errors  of  23.3%  and  60.3%  respectively. 

Figure  4.9  shows  coherence  length  trends  for  Cn*  =  10*13 
using  the  Martin  and  Flatt6  method  of  filtering.  Ten  samples 
were  obtained  for  each  data  point.  Results  were  4.07  mm  for 
the  512  x  512  array  and  4.46  mm  for  the  256  x  256  array, 
representing  errors  of  35.7%  and  48.7%  respectively. 

Figure  4.10  shows  coherence  length  trends  for  Cn*  =  10- 14 
also  using  the  Martin  and  Flatt6  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  14.94 
mm  for  the  512  x  512  array  and  20.46  mm  for  the  256  x  256 
array,  representing  errors  of  25.0%  and  71.2%  respectively. 

Figure  4.11  shows  coherence  length  trends  for  Cn*  =  10-10 
also  using  the  Martin  and  Flatt6  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  60.74 
mm  for  the  512  x  512  array  and  74.04  mm  for  the  256  x  256 
array,  representing  errors  of  27.7%  and  55.6%  respectively. 
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Figure  4.7  Calculated  Coherence  Length  Values  For 
Cn*  =  10-»«  Using  Different  Aperture  Array  Sizes 
And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.8  Calculated  Coherence  Length  Values  For 
Cn2  =  10"1B  Using  Different  Aperture  Array  Sizes 
And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.9  Calculated  Coherence  Length  Values  For 
Cn*  *  10" 1 3  Using  Different  Aperture  Array  Sizes 
And  Martin  &  Flatt6  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 


45 


40-. 


•••••  512  by  512  array 
o  o  o  o  o  256  by  256  array 
- Theoretical  coherence  length 


S  30- 

fi  : 


4> 

u 

V 

jd 

o 

u 

10- 


0  H — i — i — i — i — i — i — i — r— i — | — i — i — i — r — t — r — i — i — r — [— r~i — 1 — T — i — i — t — i — i — | 

0  100  200  300 

Aperture  Sample  Points  (pixels) 


Fl^,U»e  4,10  Calculated  Coherence  Length  Values  For 
Cn  -  10-‘«  Using  Different  Aperture  Array  Sizes 
And  Martin  &  Flatt4  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.11  Calculated  Coherence  Length  Values  For 
Cn2  =  10-18  Using  Different  Aperture  Array  Sizes 
And  Martin  &  Flatt6  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 
Representing  The  Standard  Deviation  For  Ten  Samples. 


The  most  accurate  results  were  obtained  when  the  aperture 
array  was  one-fourth  the  size  of  the  electric  field  array,  as 
seen  in  Figures  4.6  through  4.11.  For  this  reason  all 
subsequent  coherence  lengths  were  calculated  using  a  512  by 
512  array  with  a  128  by  128  aperture  array. 

Figure  4.12  shows  coherence  length  values  calculated  with 
different  values  of  turbulence  using  the  classical  method  of 
filtering.  Ten  samples  were  obtained  for  each  calculated 
value.  Although  the  calculated  coherence  length  values 
deviated  from  theoretical  values  by  approximately  30%,  a 
linear  fit  to  the  data  shows  that  the  slopes  were  correct, 
implying  a  constant  error  in  the  algorithm  due  to 
underestimating  turbulence  effects. 

Adding  a  log-normal,  random-amplitude  screen  to  the  GAUSS 
subroutine  was  attempted  to  improve  the  simulation  results  for 
coherence  lengths.  Having  both  phase  and  amplitude  screens 
present  is  equivalent  to  the  complete  Rytov  approximation. 
Equation  2.31.  Eight  samples  were  obtained  for  each 
calculated  value.  Figure  4.13  dramatically  shows  that  this 
simple  addition  produced  calculated  coherence  lengths  that 
were  within  3%  of  theoretical  values,  from  a  linear  fit  to 
the  data. 

Table  4.1  presents  coherence  lengths  calculated  by  the 
simulation  using  both  random  phase  and  amplitude  screens  in 
the  complete  Rytov  approximation.  All  data  points  were 
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Figure  4.12  Calculated  Coherence  Length  Values  As  A 
Function  Of  Cn*  With  Only  Random  Phase 
Screen  In  The  Simulation. 

Calculated  Values  Are  About  30%  High  For  All  Levels  Of 
Turbulence,  From  Linear  Fit  To  The  Data. 

Error  Bars  Represent  Standard  Deviation  For  Ten  Samples 


Calculated  coherence  length 
Theoretical  coherence  length 


Figure  4.13  Calculated  Coherence  Length  Values  As  A 
Function  of  Cn*  With  Both  Random  Phase  And  Amplitude 
Screens  In  The  Simulation. 

Calculated  Values  Are  About  3%  High  For  All  Levels  Of 
Turbulence,  From  Linear  Fit  To  The  Data. 

Error  Bars  Represent  Standard  Deviation  For  Eight  Samples 


obtained  using  a  512  by  512  array  and  128  by  128  aperture 
array. 


TABLE  4.1  COHERENCE  LENGTHS 

This  table  shows  theoretical  values  of  coherence  lengths 
calculated  from  Equation  2.39  for  eight  different  values  of 
Cn2  . 


Cn2 

theoretical  p o  (mm) 

calculated  p o  (mm) 

1  X 

10-is 

47.57 

50 

30 

2  x 

IQ-18 

31.39 

30 

20 

5  x 

lO-io 

18.11 

14 

* 

2 

1  X 

10“14 

11.95 

9 

* 

3 

2  x 

lO"14 

7.88 

7 

* 

2 

5  x 

io-*4 

4.55 

4 

1 

1  X 

10-13 

3.00 

3 

1 

2  x 

10-13 

1.98 

2.0 

± 

0.4 

E.  INTENSITY  VARIANCE  SATURATION 

The  last  checkpoint  available  for  verifying  proper 
operation  of  the  simulation  was  whether  the  code  led  to  the 
saturation  of  intensity  variance  for  increasing  turbulence. 
This  should  happen  because  the  Rytov  assumption  of  an 
exponential  random  term,  exp(^),  has  a  magnitude  bounded  by 
plus  or  minus  one.  Figure  4.14  shows  that  saturation  does 
occur  around  Cn2  =  10-14  at  a  value  of  unity,  as  expected.  The 
normalized  variance  >  1  present  at  this  level  of  turbulence 
is  due  to  a  "dancing"  beam  centroid  caused  by  low  frequency 
tilt  of  the  electric  field  as  discussed  previously.  A  decline 
to  unity  normalized  variance  appears  in  the  calculations  for 
higher  turbulence  values.  [Ref.  12] 
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Normalized 
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Figure  4.14  Normalized  Intensity  Variance  As 
A  Function  Of  Cn* . 

Saturation  Occurs  Around  Cn1  =  10*14. 

Large  Variances  Near  Saturation  Are  Due  To  Beam 
Centroid  Displacement  Cause  By  Low  Frequency  Tilt. 
Error  Bars  Represent  Standard  Deviation  For  Ten  Samples 


V.  CONCLUSIONS 


The  program  coding  discussed  in  this  thesis  simulated  the 
propagation  of  optical  waves  through  a  random  media  using  two- 
dimensional  fast  Fourier  transform  (FFT)  techniques.  The 
coding  used  Fraunhofer  propagation  throughout .  An  extension 
to  include  multi-step  Fresnel  propagation  is  straightforward. 

Aliasing  problems  in  the  FFT  were  avoided  by  taking  sample 
points  closer  together  than  the  theoretical  coherence  length 
for  a  given  turbulence  structure  parameter,  Cn2.  Tilt 
problems  in  the  FFT  were  suppressed  by  choosing  the  aperture 
array  no  larger  than  one-half  the  size  of  the  FFT  array,  with 
one-fourth  size  giving  the  best  results. 

Accuracy  of  the  simulation  was  verified  at  various 
checkpoints  within  the  coding.  Aperture  diffraction  patterns 
and  autocorrelations  were  compared  to  expected  results. 
Turbulence  effected  diffraction  patterns  and  atmospheric 
mutual  coherence  functions  were  presented  to  show  trends  for 
increasing  turbulence.  Calculated  coherence  length  values 
were  approximately  30%  larger  than  theoretical  coherence 
values  when  using  a  512  by  512  FFT  array  with  a  128  by  128 
aperture  array.  These  results  were  obtained  with  only  a 
Gaussian  distributed  random  phase  screen  in  the  simulation. 

The  addition  of  a  Gaussi?  distributed  random  amplitude 
screen  (for  the  complete  Rytov  approximation)  to  the 
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simulation  brought  calculated  coherence  values  to  within  3% 
of  theoretical  values.  This  showed  that  a  numerical  error  was 
not  present  in  the  coding,  but  that  turbulence  effects  were 
underestimated  without  the  random  amplitude  term  included. 
A  single  step,  Fraunhofer  algorithm  must  include  the  amplitude 
term  to  simulate  turbulence  effects  on  an  optical  wave 
correctly.  Unity  saturation  of  the  normalized  intensity 
variance  occurred  for  C«*  values  of  approximately  10_l«. 

Coherence  length  values  used  two  types  of  filtering.  The 
classical  method  and  the  Martin  and  Flatt§  method  gave 
comparable  results.  The  computational  efficiency  of  the 
Martin  and  Flatt6  is  significant. 

Future  endeavors  with  this  simulation  should  include  the 
incorporation  of  Fresnel  propagation  into  the  coding.  Also, 
further  studies  with  the  complete  Rytov  approximation  in  the 
simulation  code  should  be  pursued  to  obtain  better  statistical 
results . 
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APPENDIX  A 


GENERATOR  SUBROUTINE 


c 

c 

c 

c 

c 

c 

c 

c 


c 


This  is  the  random  number  generator  algoritha  from  Dr. 
Harrison's  notes,  (Ref.  7). 


Variables: 

iran... input  seed  value  (5  digit  integer) 

r . returned  randoa  number  uniformly  distributed  from  0  to  1 


iran*iran*99947 

r=0.5  +  real(iran)*2.328306e-10 


return 

end 
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APPENDIX  B.  GAUSSIAN  DISTRIBUTION  SUBROUTINE 


c - c 

subroutine  GAUSS (narray, iran) 
c 

c  This  subroutine  changes  unifornly  distributed  random  numbers, 
c  provided  by  subroutine  RAN,  to  Gaussian  distributed  randoa 
c  numbers  with  unity  variance.  This  subroutine  is  found  in 
c  Knuth,  [Ref.  8]. 
c 

c  Variables: 


c  Phaser . real  phase  screen  2-D  array 

c  phasei . imaginary  phase  screen  2-D  array 

c  narray . dimension  of  the  2-D  arrays 

c  iran . input  seed  value  for  RAN  (5  digit  integer) 

c  r . uniformly  distributed  random  numbers 


c  xl  &  x2. .. .Gaussian  distributed  random  numbers 
c 

common  /blk2/phaser (512,512) , phasei (512,512) 
c 

do  10  i=l, narray 
do  10  j=l, narray 
20  call  RAN(iran,r) 
vl=2.*r-l. 
call  RAN(iran,r) 
v2=2.*r-l. 
s=vl*vl  +  v2*v2 
if  (s.ge.1.0)  goto  20 
scale*sqrt (-2.*alog(s)/s) 
xl=vl*scale 
x2»v2*scale 
Phaser (i,j)*xl 
phasei (i, j)=x2 
10  continue 
c 

return 

end 

c - c 
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APPDDIX  C. 


SIMULATION  CODE 


c  This  code  simulates  the  propagation  of  a  monochromatic  optical 
c  wave  in  a  turbulent  atnosphere.  Two-dinensional  FFT  routines 
c  are  used  extensively  for  calculations.  Only  Fraunhofer 
c  propagation  is  present  in  the  coding.  Filtering  is  accomplished 
c  by  a  choice  of  two  nethods:  (1)  classical  and  (2)  the  Martin 
c  and  Flattd  nethod. 
c 

c  Subroutines: 


c  IMIT . Sets  electric  field  arrays  to  zero.  A  value  of  0 

c  in  the  call  initializes  real  and  imaginary  arrays, 

c  A  value  of  1  initializes  only  the  imaginary  array, 

c  SQUARE. . .Establishes  the  planar  electric  field  for  the  case 
c  of  a  square  aperture. 

c  CIRCLE. . .Establishes  the  planar  electric  field  for  the  case 
c  of  a  circular  aperture. 

c  PLOT . Gives  a  screen  plot  of  the  array  chosen  in  the  call 

c  using  EGA  graphics.  Different  colors  are  chosen  to 

c  represent  different  orders  of  magnitude  values.  Most 

c  calls  within  this  subroutine  are  MDP  Fortran-386 

c  specific. 

c  XFORM. .. .Takes  the  2-D  FFT  of  the  two  arrays  given  in  the  call, 
c  Returns  values  in  the  same  arrays.  Makes  calls  to 

c  subroutine  FFT  (1-D  FFT  routine)  several  times.  A 

c  value  of  -1  in  the  call  is  for  a  direct  transform, 

c  A  value  of  +1  in  the  call  is  for  an  inverse  transform. 

c  FFT . 1-D  FFT  routine  supplied  by  Dr.  Don  Valters.  This 

c  routine  is  used  by  subroutine  XFORM  multiple  times 

c  to  accomplish  the  2-D  transform, 

c  This  is  Dr.  Walters  FFT 

c  MAG . Calculates  the  electric  field  magnitude  and  the 

c  intensity  field  from  the  real  and  imaginary  electric 

c  field  arrays. 

c  MCFPLOT. .Calculates  the  Mutual  Coherence  Function  (MCF)  from 
c  the  intensity  field  and  then  plots  it  on  the  screen 

c  using  EGA  graphics.  Many  calls  within  this 

c  subroutine  are  MDP  Fortran-386  specific, 

c  GAUSS. .. .Creates  the  phase  screens  from  Gaussian  distributed 
c  random  numbers  with  unit  variance.  Makes  calls 

c  to  subroutine  RAN  for  uniformly  distributed  random 

c  numbers  then  transforms  them  to  a  Gaussian  distribution 

c  Elements  of  this  subroutine  provided  by  Knuth.  [Ref.  8] 

c  RAN . Generates  uniformly  distributed  random  numbers.  This 

c  algorithm  from  Harrison  [Ref.  7] 
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c  FILTER... Filters  the  phase  screens  according  to  the  Kolmogorov 

c  power  law.  A  shuffling  of  array  values  is  necessary 

c  for  proper  filtering. 

c  MESH . Meshes  the  random  phase  screen  with  the  electric 

c  field, 

c 

c  Variables: 

c  fieldr. ..Real  2-D  electric  field  array, 
c  f ieldi. . .Imaginary  2-D  electric  field  array, 

c  fieldro. .Real  2-D  planar  electric  field  at  aperture, 
c  phaser. . .Real  2-D  phase  screen  array, 
c  phasei ... Imaginary  2-D  phase  screen  array. 

c  field*... 2-D  array  that  represents  the  electric  field  magnitude, 

c  fieldm2..2-D  array  that  represents  the  intensity  field, 

c  narray. . .Dimension  of  electric  field  array 

c  nfield. . .Dimension  of  aperture  array 

c  iran . Input  seed  for  random  number  generator  (5  digit  integer). 

c  cn2 . Refractive  index  structure  parameter. 

c  wvl . Wavelength  of  monochromatic  electromagnetic  wave. 

c  delz . Physical  distance  between  aperture  plane  and  image  plane. 

c  width _ Physical  width  of  aperture  plane. 

c  delx . Physical  distance  between  sample  points  in  aperture 

array. 

c  rhonot .. .Coherence  length  calculated  by  interpolation  between 
c  points  of  MCF  curve.  « 

c  re . Real  1-D  array  used  by  subroutine  FFT. 

c  rim . Imaginary  1-D  array  used  by  subroutine  FFT. 

c  fphaser. .Real  2-D  array  used  in  subroutine  FILTER  that  represents 
c  the  phase  screen  in  a  more  convenient  fora  for 

c  filtering  and  viewing. 

c  fphasei. .Imaginary  2-D  array  used  in  subroutine  FILTER  that 
c  represents  the  phase  screen  in  a  more  convenient  for* 

c  for  filtering  and  viewing, 

c 
c 

common  /blkl/fieldr (512,512) ,f ieldro(512,512) , field! (512, 512) 
common  /blk2/phaser (512,512) ,phasei (512,512) 
common  /blk3/fielda(512,512) ,field*2 (512,512) 
c 

c  Initialize  arrays 
c 

1  call  WIT  (narray,  0) 
c 

c  Input  section 
c 

write(*,*)  'Enter  dimension  of  array  that  required.' 
read(*,*)  narray 

write(*,*)  'Enter  dimension  of  planar  electric  field.’ 
write(*,*)  '(Pixel  width  of  the  aperture)' 
read(*,*)  nfield 

2  write (*,*)  'Choose  aperture  shape:  1)  SQUARE' 

write(*,*)  '  2)  CIRCLE' 
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c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


read(*,*)  i choice 

if  (ichoice  .It.  1  .or.  ichoice  .gt.  2)  then 
write(*,*)  'Try  again! * 
goto  2 
endif 

write(*,*)  'Enter  the  seed  for  the  random  number  generator.' 
write(*,*)  '(Value  must  be  a  five  digit  integer)' 
read(*,*)  iran 

write(*,*)  'Enter  the  value  of  Cn  squared.' 
read(*,*)  cn2 

write(*,*)  'Enter  the  wavelength  of  light.' 
read(*,*)  wvl 

write(*,*)  'Enter  the  distance  from  aperture  to  observer.' 

read(*,*)  delz 

delz=1000. 

write(*,*)  'Enter  the  width  of  the  aperture  in  meters.' 
read(*,*)  width 

delx=width/real (nf ield) 


Load  arrays  with  desired  planar  electric  field 

if  (ichoice  .eq.  1)  call  SQUARE (narray,nf ield) 
if  (ichoice  .eq.  2)  call  CIRCLE (narray,nf ield) 

Plot  planar  electric  field 

call  PLOT(fieldr,narray,l) 

Take  Fourier  transform  of  the  field 


call  XFORM (f ieldr , f ieldi , narr ay , delx , -1 . ) 

Calculate  magnitude  of  transformed  field  then  plot  it 
call  KAG(narray) 
call  PLOT ( field*, narr ay, 1) 

Set  imaginary  portion  of  field  to  zero 
call  INIT(narray,l) 

Take  Fourier  transform  of  field  intensity 
call  XF0RM(fieldm2, f ieldi, narr ay, delx, +1 . ) 

Calculate  and  plot  MCF  of  aperture 
call  MCFPLOT (narray, delx, nf ield, 0) 


Load  phase  screen  arrays  with  c  ssian  random  numbers 
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c 

call  GAUSS (narray, iran) 
c 

c  Transfora  phase  screen  to  frequency  space 
c 

call  XFORM (phaser , phasei , narray , sqrt (real (narray) ) , -1 . ) 
c  call  PLOT (phaser, nar ray ,1) 

c  call  PLOT (phasei, narray, 1) 

c 

c  Filter  phase  screen  using  Kolmogorov  spectrum  idea 
c 

call  FILTER (narray, cn2,wvl,delx,delz) 
c  call  PLOT (phaser, narray, 1) 

c  call  PLOT (phasei, nar r ay, 1) 

c 

c  Transform  phase  screen  back  to  real  space 
c 

call  XFORN (phaser, phasei, narray, sqrt (real (narray )*delx) ,+l.) 
c  call  PLOT (phaser, nar r ay, 1) 

c  call  PLOT (phasei, nar ray, 1) 

c 

c  Mesh  phase  screen  with  planar  electric  field 
c 

call  IlflT (narray,  1) 
call  MESH (narray) 
c 

call  PLOT(f ieldr, narray, 1) 
c 

c  Take  fourier  transform  of  this  electric  field 
c 

call  XFORM(f ieldr ,f ieldi, narray, delx,-l. ) 
c 

call  MAG (narray) 

c 

call  PLOT(f ieldm, narray, 1) 
c 

c  Reset  imaginary  portion  of  field  to  zero 
c 

call  INIT (narray, 1) 
c 

c  Take  Fourier  transform  to  get  diffraction  pattern  w/  turbulence, 
c 

call  XFORM ( f ieldm2 , f ieldi , narray , delx , +1 . ) 
c 

c  Calculate  atmospheric  MCF  and  plot  it. 
c 

call  MCFPLOT (narray, delx, nfield,l) 
c 

goto  1 

stop 

end 


6' 


c 


subroutine  INIT(narray,k) 


« 


c 


c 


c 

c 


c 


c 


c 

c 

c 


c 


c 


c 


c 


couod  /blkl/f ieldr (512 f  512) ,  fieldro (512 , 512) , f ieldi  (512,512) 
do  10  i=l,narray 
do  10  j*l,narray 
if  (k  .eq.  0)  then 
fieldr (i, j)=0.0 
fieldro(i, j)=0.0 
end  if 

fieldi(i, j)=0.0 
10  continue 

return 

end 

- c 

subroutine  MAG(narray) 

coaaon  /blkl/f ieldr (512,512) , fieldro (512, 512) ,fieldi(512,512) 
coBBon  /blk3 /fields (51 2, 512) ,field»2(512,512) 

do  10  i=l,narray 
do  10  j=l,narray 

fields2(i, j)=fieldr(i, j)**2  +  fieldi(i, j)**2 
fieldm(i, j)=sqrt (fieldn2(i, j) ) 

10  continue 

return 

end 

- c 

subroutine  SQUARE (narray,nf ield) 

cobbod  /blkl/f ieldr (512,512) , fieldro (512, 512) ,f ieldi (512,512) 

do  10  i=narray/2+l-nfield/2,narray/2+nf ield/2 
do  10  j®narray/2+l-nf ield/2, narray/2+nf ield/2 
fieldr (i, j)=1.0 
fieldro(i, j)*1.0 
10  continue 

return 

end 

- c 

subroutine  CIRCLE(narray ,nf ield) 

cowaon  /blkl/f ieldr (512, 512) ,fieldro(512,512) , f ieldi (512,512) 

narray2*narray/2 

nfield2*nf ield/2 

do  10  i=narray2+l-nfield2,narray2+nfield2 
do  20  j=narray2+l-nf ield2,narray2+nfield2 
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radius2* (real (i) -real (narray2) ) **2  + 

*  (real (j) -real (narray2))**2 

if  (radius2  .It.  real(nfield2)**2)  then 
fieldr (i, j)=1.0 
fieldroti, j)*1.0 
endif 

20  continue 
10  continue 
c 

return 

end 

c - c 

subroutine  PLOT(field,narray,k) 
c 

dimension  field (512, 512) 
dimension  ndex(20) 

data  ndex/4, 4, 12, 12, 14, 14. 10, 10, 2, 2, 3, 3, 9, 9, 1,1, 8, 8, 0,0/ 
c 

f max3 0.0 

do  100  i=l,narray 
do  100  j«l,narray 
x=field(i, j) 
if  (x.gt.fmax)  then 
fmax^x 
endif 

100  continue 
c 

call  set_yideo_node(16) 
call  eaa_set_mode_4 
c 

do  110  i=l,narray/k 
do  110  j®l,narray/k 
ix=i 
iy»j 

xmax=alogl0 (f ield (i , j ) /f max) 
index»8*abs (xmax) +1 
if  ( index. ge. 21)  then 
icolor*0 
else 

icolor=ndex( index) 
endif 

call  ega_put_pixel(ix,iy,icolor) 

110  continue 
c 

call  pause 

call  set_video_mode (3) 
c 

return 

end 

c - c 

subroutine  MCFPLOT(narray,delx,nf ield.iaprture) 
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common  /blk3/f ieldm(512,512)  ,fieldm2(512,512) 
dimension  apmcf (513) ,ymcf (513) 
logical  flag 
efold=0. 36787944 
rhonot=0. 
flag*. true, 
c 

do  10  i*l,513 

if  (iaprture  .eq.  0)  then 
apmcf (i)=0.0 
endif 

ymcf (i)=0.0 
10  continue 
c 

f»cf=0.0 

do  20  i=l,narray 
do  20  j*l,narray 
xmcf=f ieldm2(i, j) 
if  (xmcf .gt.fmcf )  then 
fmcf=xmcf 
endif 

20  continue 
c 

do  30  i=l,nfield 

if  (iaprture  .eq.  0)  then 
apmcf (i) =f ield»2 (1 ,  i) It mcf 
ymcf (i)=apmcf (i) 
else 

ymcf (i)=field»2(l,i) /(apmcf (i)*fncf ) 
endif 

30  continue 
c 

do  40  i=l,nfield 

if  (iaprture. eq.l)  then 
if  (flag)  then 

if  (ymcf (i) .lt.efold)  then 
rhonot=real (i-1) *delx 

*  +  (efold-ymcf (i) )*delx/ (ymcf (i-1) -ymcf (i) ) 

flag*. false, 
endif 
endif 
endif 

40  continue 
c 

call  set_video_mode(16) 
call  ega_set_node_4 
call  locate(20,l) 
if  (iaprture  .eq.  0)  then 

call  write_string( 'Aperture  MCF  vs.  Length') 
else 

call  write_string{ 'Atmospheric  MCF  vs.  Length') 
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endif 

call  locate (1,1) 
call  write_string( '1.0' ) 
call  locate (25, 24) 
call  write_string( 'length' ) 
c  call  locate (37 ,24) 

c  call  write_atring(nfield) 

call  ega_draw_line (40 ,20,40,320,15) 
call  ega_draw_line (40,320,600,320,15) 
c 

ix*40 

iy»20 

do  50  isl,nfield 
ix2=40+560*i/nfield 
iy2*20+int (300. * (1 . -ymcf (i+1) ) ) 
call  ega_dr aw_line ( ix , iy , ix2 , iy 2 , 15 ) 
ix-ix2 
iy=iy2 
50  continue 
c 

call  pause 

call  set_videojaode{3) 
c 

write(*,*)  'The  coherence  length  is  ',rhonot,'  neters.’ 
c 

return 

end 


subrout ine  XFORM ( £ ieldr , f ieldi , nar r ay , £  f  t nor m , s ign ) 

connon  /blk4/re(512) ,ria(512) 
dimension  f ieldr (512,512) ,f ieldi (512,512) 
data  re/512*0./,rim/512*0./ 

«=int (alog (real (narray) ) /alog (2 . ) ) 

do  30  i=l, narray 
do  40  j=l,narray 
re(j)*fieldr (i, j) 
rim(j)*£ieldi(i, j) 

40  continue 

call  FFT (a, f f tnora, sign) 
do  50  j=l, narray 
fieldr (i, j)=re(j) 
fieldi(i, j)=rin(j) 

50  continue 
30  continue 

do  60  j=l, narray 
do  70  i=l, narray 
re(i)*fieldr (i, j) 
rin(i)*£ieldi(i, j) 
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70  continue 

call  FFT(s,fftnor«,sign) 
do  80  i*l,narray 
fieldr (i, j)*re(i) 
fieldi{i, j)«ri»(i) 

80  continue 
60  continue 

return 

end 


subroutine  FFT(a,£ftnora,sign) 

couon  /blk4/re(512)  ,ris(512) 

pi=3.141592653589792*sign 

n=2**« 

nl=n-l 

3=1 

do  200  i=l,nl 
if  (i.lt.j)  then 
t=re(j) 
re(j)=re(i) 
re(i)=t 
t=ria(j) 
ria(j)=ria(i) 
ria(i)=t 
endif 
k=n/2 

do  201  while  (k.lt.j) 

3=3~k 

k=k/2 

201  continue 
j=j+k 

200  continue 
le=l 

do  202  1~1 f n 
lel=le 
le=le+le 
ure=l. 
uia=0 . 
angspi/lel 
wre*cos(ang) 
wia=sin(ang) 
do  203  j=l,lel 
do  204  i=j,n,le 
ip*i+lel 

tre=re(ip) *ure-ria(ip) *ui« 

tia=re (ip) *uia+ria (ip) *ure 

re(ip)=re(i)-tre 

ria(ip)«ria(i)-tia 

re(i)=re(i)+tre 

ria(i)*ria(i)+tia 
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204  continue 
t*ure*wre-uin*wi« 
ui«*ure*wi»+ui**wre 
ure=t 

203  continue 
202  continue 

if  (sign. gt. 0.0)  then 
pts=1.0/real(n*ff tnorn) 
else 

pts=f ftnorn 
endif 

do  205  i«l,n 
re(i)=re(i) *pts 
ri«(i)=ri«(i) *pts 

205  continue 

return 

end 


subroutine  GAUSS (narray,iran) 

couon  /blk2/phaser (512,512)  ,phasei (512,512) 

do  10  i=l,narray 
do  10  j=l,narray 
20  call  RAN(iran.r) 
vl-2.*r-l. 
call  RAN(iran,r) 
v2*2.*r-l. 
s=vl*vl  +  v2*v2 
if  (s.ge.1.0)  goto  20 
scale=sqrt (-2.*alog(s)/s) 
xl=vl*scale 
x2=v2*scale 
Phaser (i, j )*xl 
phasei(i, j)=x2 
10  continue 

return 

end 


subroutine  RAN(iran,r) 

kran*99947 

iran*iran*kran 

r*0.5  +  real(iran)*2.328306e-10 

return 

end 


subroutine  FILTER (nar ray , cn2 , wvl , delx, delz) 


c 


c 


c 
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coaaon  /blk2/phaser (512,512) ,phasei(512,512) 

dimension  f phaser (512,512) ,fphasei(512,512) 

dimension  aaa(512,512) 

pi*3. 141592653589792 

tpi*2.*pi 

narray2=narray/2 

npivot=narray2+l 

power=-ll./6. 

factor=sqrt (1.303* (tpi**3)*cn2*delz)/wvl 
f actor2= ( tpi/real (delx*narray) ) **power 

do  10  iel,narray 
do  10  j»l,narray 
f Phaser (i, j)=0. 
fphasei (i, j)*0. 

10  continue 

do  20  i=l,narray 
do  20  j=l,narray 

if  (i.le.npivot)  then 
if  (j.le.npivot)  then 

f Phaser (i-l+narray2 , j-l+narray2) =phaser (i , j ) 
fphasei (i-l+narray2,j-l+narray2)=phasei (i, j) 
else 

f Phaser (i-l+narray2 , j-l-narray2) =phaser (i , j ) 
fphasei (i-l+narray2 , j-l~narray2) *phasei (i , j ) 
endif 
else 

if  (j.le.npivot)  then 

f phaser ( i-l-narray2 , j -l+narray2) ®phaser ( i , j ) 
fphasei (i-l-narray2, j-l+narray2)=phasei (i, j) 
else 

f phaser (i-l-narray2 , j-l-narray2) =phaser (i , j ) 
fphasei (i-l-narray2, j-l-narray2)=phasei (i, j) 
endif 
endif 
20  continue 

do  30  i=l,narray 
do  30  j=l,narray 

freq=sqrt (real (i-narray2) **2  +  real (j-narray2) **2) 
if  (i.eq.narray2.and. j .eq.narray2)  then 
f Phaser (i, j)=0. 
fphasei (i, j)*0. 
else 

f Phaser (i, j) *f phaser (i, j) *f actor *f actor 2* (f req) **power 
f phasei (i , j) =f phasei (i , j) *f actor*f actor 2* (f req) **pover 
endif 

30  continue 
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do  40  i«l,narray 
do  40  j>l,narray 

if  (i.lt.narray2)  then 
if  (j .lt.narray2)  then 
Phaser ( i+l+narr ay2 , j+l+narray2) =f phaser (i ,  j ) 
phasei ( i+l+narray 2 , j+l+nar ray2 ) *f phasei ( i , j ) 
else 

Phaser (i+l+narray2, j+l-narray2) =f phaser (i , j) 
phasei ( i+l+narr ay2 , j+l-narray2) *f phasei ( i , j ) 
endif 
else 

if  (j . It.narray2)  then 
Phaser (i+l-narray2 , j+l+narray 2) «f phaser (i , j ) 
phasei ( i+l-narray2 , j+l+narray2) =f phasei ( i , j ) 
else 

Phaser (i+l-narray2 , j+l-narray2) =f phaser (i , j ) 
phasei (i+l-narray2, j+l-narray2)sf phasei (i,  j) 
endif 
endif 

40  continue 
c 

return 

end 

c - c 

subroutine  NESH(narray) 
c 

couon  /blkl/f ieldr  (512 , 512) , fieldro(512, 512)  ,fieldi (512,512) 
conon  /blk2/phaser  (512,512) , phasei (512, 512) 

c 

do  10  i=l,narray 
do  10  j=l,narray 

ercosphi=f ieldro ( i , j ) *cos (phaser { i , j ) ) *exp (phasei ( i , j ) ) 
eicosphi*f ieldi (i, j)*cos (phaser (i, j) )*exp (phasei (i, j) ) 
ersinphi=f ieldro (i , j ) *sin (phaser (i , j ) ) *exp (phasei (i , j ) ) 
eisinphi=f ieldi (i,j)*sin (phaser (i, j) ) *exp (phasei (i, j) ) 
f ieldr (i , j ) =ercosphi-eisinphi 
f ieldi (i,j)*ersinphi+eicosphi 
10  continue 
c 

return 

end 
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